/*rtnode.cpp
  this file implements class RTNode*/

#include <math.h>
#include <stdlib.h>
#include <string.h>
#include "rtnode.h"
#include "rtree.h"
#include "entry.h"
#include "../blockfile/blk_file.h"
#include "../blockfile/cache.h"
#include "../linlist/linlist.h"
//------------------------------------------------------------
RTNode::RTNode(RTree *rt)
  //use this ructor to create a new node on disk.
{
    char *b;
    int header_size;
    Entry * d;
    int i;

    my_tree = rt;
    dimension = rt->dimension;
    num_entries = 0;
	dirty = TRUE;

    d = new Entry();
	d -> init_entry(dimension, NULL);
    header_size = sizeof(char) + sizeof(int);  // level + num_entries
    capacity = (rt -> file -> get_blocklength() - header_size) / d -> get_size();
    delete d;

    entries = new Entry[capacity];
    for (i = 0; i < capacity; i++)
		entries[i].init_entry(dimension, rt);

	//assign a new block on the disk
    b = new char[rt -> file -> get_blocklength()];
    block = rt -> file -> append_block(b);  //add a new block to the index file
    delete [] b;
}
//------------------------------------------------------------
RTNode::RTNode(RTree *rt, int _block)
  //use this ructor to restore a node from the disk.
{
    char *b;
    int header_size;
    Entry * d;
    int i;

    my_tree = rt;
    dimension = rt->dimension;
    num_entries = 0;
	dirty = FALSE;

    d = new Entry();
	d -> init_entry(dimension, NULL);
    header_size = sizeof(char) + sizeof(int);
    capacity = (rt -> file -> get_blocklength() - header_size) / d -> get_size();
    delete d;

    entries = new Entry[capacity];
    for (i = 0; i < capacity; i++)
		entries[i].init_entry(dimension, rt);
    
    block = _block;
    b = new char[rt -> file -> get_blocklength()];
    if (rt -> cache == NULL) // no cache
        rt -> file -> read_block(b, block);
    else
        rt -> cache -> read_block(b, block, rt);

    read_from_buffer(b);
    delete [] b;
}
//------------------------------------------------------------
RTNode::~RTNode()
{
    char *b;

    if (dirty)
    {
        b = new char[my_tree->file->get_blocklength()];
        write_to_buffer(b);

        if (my_tree->cache == NULL) // no cache
            my_tree->file->write_block(b, block);
        else
            my_tree->cache->write_block(b, block, my_tree);

        delete [] b;
    }

    delete [] entries;
}

/*****************************************************************
this function returns the best sub-tree to insert the new entry

para:
newe: the entry to be inserted

*****************************************************************/

int RTNode::choose_subtree(Entry *_newe)
{
	Entry *hypo_e=new Entry(dimension, my_tree); //the hypothetical entry
    float fmin = MAXREAL;
    float amin = MAXREAL;

	int minindex=-1;
    for (int i = 0; i < num_entries; i++)
    {
		*hypo_e=entries[i];
		future_mbr(hypo_e->bounces, hypo_e->velocity, my_tree->time, dimension);
		float a=IntgArea(hypo_e->bounces, hypo_e->velocity, 0, my_tree->T);
		hypo_e->enlarge(hypo_e, _newe);

        // calculate area and area enlargement
        float f=IntgArea(hypo_e->bounces, hypo_e->velocity, 0, my_tree->T)-a;
			//Explanation: note hypo_e stores the extent at insertion time

    	if ((f<fmin) || (f==fmin && a<amin))
    	{
    	    minindex = i;
        	fmin = f;
        	amin = a;
        }
    }
	delete hypo_e;
    return minindex;
}

//------------------------------------------------------------
void RTNode::enter(Entry *de)
  //note that de will be deleted after being entered.
{
    if (num_entries > (capacity-1))
        error("RTNode::enter: called, but node is full", TRUE);

	past_mbr(de->bounces, de->velocity, my_tree->time, dimension);
    entries[num_entries] = *de;

    num_entries++;

	dirty = true;

    de->son_ptr = NULL;
    delete de;
}
//------------------------------------------------------------
bool RTNode::FindLeaf(Entry *_q)
{
	bool ret=false;
	if (level==0)
	{
		for (int i=0; i<num_entries; i++)
		{
			if (_q->son==entries[i].son)
			{
				printf("\nfind the record in node %d, ", block);
				return true;
			}
		}
	}
	else
	{
		for (int i = 0; i < num_entries; i++)
		{
			float mbr1[4], mbr2[4];
			memcpy(mbr1, _q->bounces, sizeof(float)*4);
			memcpy(mbr2, entries[i].bounces, sizeof(float)*4);
			future_mbr(mbr1, _q->velocity, my_tree->time, 2);
			future_mbr(mbr2, entries[i].velocity, my_tree->time, 2);
			float *intv=overlapRect(2, mbr1, mbr2);
//			float *intv=my_tree->moving_sect(&(entries[i]), _q, my_tree->time, my_tree->time);
			if (intv)
			{
				RTNode *succ = entries[i].get_son();
				bool cret=succ->FindLeaf(_q);
				if (cret)
					printf(" ->entry %d", i);
				ret=ret|cret;
				entries[i].del_son();
			}
		}
	}
	return ret;
}

/*****************************************************************
This function initiates a float array and return the mbr
of the node at the current time of the TPR-tree. this function also
tightens the mbr at the current time

*****************************************************************/

float* RTNode::get_mbr()
{
    int i, j;
    float *mbr;

    mbr = new float[2*dimension];
    for (i = 0; i < 2*dimension; i ++ )
        mbr[i] = entries[0].bounces[i]+entries[0].velocity[i]*my_tree->time;

    for (j = 1; j < num_entries; j++)
    {
    	for (i = 0; i < 2*dimension; i += 2)
    	{
    	    mbr[i]   = min(mbr[i],   entries[j].bounces[i]+entries[j].velocity[i]*my_tree->time);
    	    mbr[i+1] = max(mbr[i+1], entries[j].bounces[i+1]+entries[j].velocity[i+1]*my_tree->time);
        }
    }

    return mbr;
}

/*****************************************************************
This function initiates a float array and return the velocity mbr
of the node

*****************************************************************/

float* RTNode::get_vmbr()
{
    int i, j;
    float *vmbr;

    vmbr = new float[2*dimension];
    for (i = 0; i < 2*dimension; i ++ )
        vmbr[i] = entries[0].velocity[i];

    for (j = 1; j < num_entries; j++)
    {
    	for (i = 0; i < 2*dimension; i += 2)
    	{
    	    vmbr[i]   = min(vmbr[i],   entries[j].velocity[i]);
    	    vmbr[i+1] = max(vmbr[i+1], entries[j].velocity[i+1]);
        }
    }

    return vmbr;
}

/*****************************************************************
This function initiates a float array and return the velocity mbr
of the node considered against boundary
Wrote for new type of conservative bounding rectangle

*****************************************************************/

float* RTNode::get_vmbr(float *_mbr)
{
    int i, j, H;
    float *vmbr, position; 

	H=my_tree->T; //- my_tree->time; //0907
    vmbr = new float[2*dimension];
    for (i = 0; i < 2*dimension; i ++ )
	{
		position = entries[0].bounces[i]+entries[0].velocity[i]*my_tree->time;
		vmbr[i] = (entries[0].velocity[i]*H-_mbr[i] +position)/H;
	}

    for (j = 1; j < num_entries; j++)
    {
    	for (i = 0; i < 2*dimension; i += 2)
    	{
    	    position = entries[j].bounces[i]+entries[j].velocity[i]*my_tree->time;
			vmbr[i]   = min(vmbr[i],   (entries[j].velocity[i]*H-_mbr[i] +position)/H);
    	    position = entries[j].bounces[i+1]+entries[j].velocity[i+1]*my_tree->time;
			vmbr[i+1] = max(vmbr[i+1], (entries[j].velocity[i+1]*H-_mbr[i+1] +position)/H);
		}
    }
	
	return vmbr;
}

/*****************************************************************
Use this function to insert an entry into the node
para:
d: the entry to be inserted. Note that the MBR of d must be at the 
current time
sn: the pointer to the new node (if created from split)
*****************************************************************/

R_OVERFLOW RTNode::insert(Entry *d, RTNode **sn)
{
	R_OVERFLOW ret;
    if (level > d->level)
    {
		//get the subtree to insert---------------------------
		int follow=choose_subtree(d); 
		//insert it-------------------------------------------
		RTNode *succ=entries[follow].get_son();
		RTNode *new_succ; //point to the new node (if created from split)
		//09/21
		/* if (level==d->level+1==1)  //if this node is parent of d
		{
			if ((d->velocity[0]+d->velocity[1])*(entries[follow].velocity[0]+entries[follow].velocity[1])>0 && (d->velocity[2]+d->velocity[3])*(entries[follow].velocity[2]+entries[follow].velocity[3])>0)
				R_OVERFLOW ret=succ -> insert(d, &new_succ); // new_succ is not null if after insert, a split occured
			else
			{
				new_succ=new RTNode(my_tree);
				new_succ->level = d->level;
				if (d->level==0)
					my_tree->num_of_dnodes++;
				else
					my_tree->num_of_inodes++;
				new_succ->enter(d);
				ret = SPLIT;
			}
		}
		else
		{
		*/
			R_OVERFLOW ret=succ -> insert(d, &new_succ); // new_succ is not null if after insert, a split occured
		//}
				//adjust the parent entry-----------------------------
		update_parent_entry(follow);
		//----------------------------------------------------
		entries[follow].del_son();
		//----------------------------------------------------
		if (ret == SPLIT)  // check whether child node has been split or not
		{
			if (num_entries == capacity)
				error("RTNode::insert--maximum capacity violation", true);
			//init a new entry for the new node---------------
			entries[num_entries].son=new_succ->block;
			entries[num_entries].son_ptr=new_succ;
			update_parent_entry(num_entries);
			entries[num_entries].del_son();
			num_entries++;
			//split the node if necessary --------------------
			if (num_entries == (capacity - 1)) //to be split
			{
				*sn=new RTNode(my_tree);
				my_tree->num_of_inodes++;
				(*sn)->level=level;
				split(*sn);   
				ret = SPLIT;
			}
			else
				ret = NONE;
			//------------------------------------------------
		}
		dirty = true;
		return ret;
    }
    else 
    {
        if (num_entries == capacity)
        	error("RTDataNode::insert--maximum capacity violation", TRUE);
     // insert new entry into the leaf node found by choosesubtree
		enter(d);
		//split the node if necessary ----------------------------
		if (num_entries == (capacity - 1)) //to be split
		{
	//====================================================
					//OVERFLOW TREATMENT//	

		//REINSERT-OVERFLOW TREATMENT	

			if (level==0 && my_tree->re_level[0] == false && my_tree->root_ptr->level>0)
            {
                float *mbr=get_mbr();
				float *vmbr=get_vmbr(mbr);
				past_mbr(mbr, vmbr, my_tree->time, 2);
				
        	    TPRSortMbr *sm = new TPRSortMbr[num_entries];
        	    for (int i = 0; i < num_entries; i++)
        	    {
            		sm[i].index = i;
            		sm[i].dimension = dimension;
					sm[i].st=my_tree->time;
					sm[i].ed=my_tree->time+my_tree->T;
					sm[i].bounces=entries[i].bounces;
					sm[i].velocity=entries[i].velocity;
            		sm[i].centermbr=mbr;
            		sm[i].centervmbr=vmbr;
                }
                //sort the entries by their center distances------
                qsort(sm, num_entries, sizeof(TPRSortMbr), IntgSortCenMbr);
                // reinsert about 30%
				int last_cand = (int) ((float)num_entries * 0.30);

                // copy the nearest candidates to new array ------
				Entry *new_entries = new Entry[capacity];
			    for (int i = 0; i < capacity; i++)
					new_entries[i].init_entry(dimension, my_tree);
                for (int i = 0; i < num_entries - last_cand; i++)
    	            new_entries[i] = entries[sm[i].index];
                // insert candidates into reinsertion list -------
                for (int i=0; i<num_entries; i++)
                {
					//recalculate for current time
					future_mbr(entries[sm[i].index].bounces, entries[sm[i].index].velocity,
						my_tree->time, dimension);
					Linkable *new_link=entries[sm[i].index].gen_Linkable();
					new_link->level=0;
                    my_tree->re_data_cands->insert(new_link);
                }
				//replace entries in this node with the new ones--
				delete [] entries;
        	    entries = new_entries;
				//------------------------------------------------
				delete []vmbr; delete []mbr; delete []sm;
        	    my_tree->re_level[0] = true;
        	    num_entries -= last_cand;
        	    dirty = true;
                return REINSERT;
        	}
		//SPLIT-OVERFLOW TREATMENT
			else
			{
				*sn=new RTNode(my_tree);
				(*sn)->level = level;
				my_tree->num_of_dnodes++;
				split(*sn); 
   				ret=SPLIT;
			}
			//END OF OVERLOW TREATMENT//
	//=================================================
        }
        else
            ret=NONE;
		dirty=true;
		return ret;
    }
}
/*****************************************************************
this function answers a window query recursively

para:
q: the moving query
st: the starting time of the itnerval
ed: the ending time
res: a linked list storing the result

*****************************************************************/

void RTNode::rangeQuery(Entry *_q, float _st, float _ed, int& _rsltcnt)
//void RTNode::rangeQuery(Entry *_q, float _st, float _ed, SortedLinList *res)
{
    for (int i = 0; i < num_entries; i++)
    {
//if (entries[i].son==8597)
//printf("testing...\n");
        float *intv=my_tree->moving_sect(&(entries[i]), _q, _st, _ed);
		if (intv)
        {
			delete []intv;
            if (level == 0)
            {
//enable this line if you want to see the actual objects----------
//              Linkable *copy;
//				copy = entries[i].gen_Linkable();
//				res -> insert(copy);  
//----------------------------------------------------------------
				_rsltcnt++;
//testing---
//printf("id=%d, (%f, %f), velocities=(%f, %f)\n", entries[i].son, entries[i].bounces[0], entries[i].bounces[2],
//	   entries[i].velocity[0], entries[i].velocity[2]);
            }
            else
            {
                RTNode *succ = entries[i].get_son();
                //succ -> rangeQuery(_q, _st, _ed, res);
				succ -> rangeQuery(_q, _st, _ed, _rsltcnt);
				entries[i].del_son();
            }
        }
    }
}
//------------------------------------------------------------
void RTNode::read_from_buffer(char *buffer)
{
    int i, j, s;

    // Level
    memcpy(&level, buffer, sizeof(char));
    j = sizeof(char);

    // num_entries
    memcpy(&num_entries, &(buffer[j]), sizeof(int));
    j += sizeof(int);

    s = entries[0].get_size();
    for (i = 0; i < num_entries; i++)
    {
    	entries[i].read_from_buffer(&buffer[j]);
    	j += s;
    }
}


///*
int RTNode::split(float **_mbr, float **_vmbr, int **_distribution)
{
	int k,s,l,i;
    int n = num_entries;
    int m1 = (int) ceil((float)n * 0.40);  //minimum number of entries in a node

    TPRSortMbr *sml = new TPRSortMbr[n];
    TPRSortMbr *smu = new TPRSortMbr[n];
	float *rxmbr = new float[2*dimension];  //the mbr of the first node
//	float *rymbr = new float[2*dimension];  //the mbr of the second node
	float *rxvmbr=new float[2*dimension];  //the vmbr of the first node
//	float *ryvmbr=new float[2*dimension];  //the vmbr of the second node

    // choose split axis------------------------------------------
    float minmarg = MAXREAL; 
	int split_axis=-1;
    for (int i = 0; i < dimension; i++)    // for each axis
    {
		//first check spatial dimensions #########################
		//assembly sml and smu for sorting------------------------
        for (int j = 0; j < n; j++)
        {
            sml[j].index=smu[j].index=j;
            sml[j].dimension=smu[j].dimension=i;
            sml[j].bounces=smu[j].bounces=_mbr[j];
			sml[j].velocity=smu[j].velocity=_vmbr[j];
			sml[j].st=smu[j].st=my_tree->time;
			sml[j].ed=smu[j].ed=my_tree->time;
        }
        // Sort by lower and upper -------------------------------
      	qsort(sml, n, sizeof(TPRSortMbr), IntgSortHighMbr);
        qsort(smu, n, sizeof(TPRSortMbr), IntgSortLowMbr);
		//--------------------------------------------------------
        float marg = 0.0;
        // for all possible distributions of sml------------------
        for (int k = 0; k < n - 2 * m1 + 1; k++)
        {
			//get the extents of the first node-------------------
			for (int s = 0; s < 2 * dimension; s += 2)
			{
				rxmbr[s] =    MAXREAL;  rxvmbr[s]=MAXREAL;
				rxmbr[s+1] = -MAXREAL;  rxvmbr[s+1]=-MAXREAL;
			}
            for (int l = 0; l < m1 + k; l++)
            {
				for (int s = 0; s < 2*dimension; s += 2)
				{
					rxmbr[s] =   min(rxmbr[s],   sml[l].bounces[s]+sml[l].velocity[s]*my_tree->time);
					rxmbr[s+1] = max(rxmbr[s+1], sml[l].bounces[s+1]+sml[l].velocity[s+1]*my_tree->time);
					rxvmbr[s] =   min(rxvmbr[s],   sml[l].velocity[s]);
					rxvmbr[s+1] = max(rxvmbr[s+1], sml[l].velocity[s+1]);
				}
			}
			//add the margin--------------------------------------
			marg += IntgMargin(rxmbr, rxvmbr, 0, my_tree->T);
			//extents of the second node--------------------------
			for (int s = 0; s < 2 * dimension; s += 2)
			{
				rxmbr[s] =    MAXREAL; rxvmbr[s]=MAXREAL;
				rxmbr[s+1] = -MAXREAL; rxvmbr[s+1]=-MAXREAL;
			}
            for (int l = 0; l < n; l++)
            {
				for (int s = 0; s < 2 * dimension; s += 2)
				{
					rxmbr[s] =   min(rxmbr[s],   sml[l].bounces[s]+sml[l].velocity[s]*my_tree->time);
					rxmbr[s+1] = max(rxmbr[s+1], sml[l].bounces[s+1]+sml[l].velocity[s+1]*my_tree->time);
					rxvmbr[s] =   min(rxvmbr[s],   sml[l].velocity[s]);
					rxvmbr[s+1] = max(rxvmbr[s+1], sml[l].velocity[s+1]);
				}
            }
			//add the margin--------------------------------------
			marg += IntgMargin(rxmbr, rxvmbr, 0, my_tree->T);
			//----------------------------------------------------
        }

        // for all possible distributions of smu------------------
       	for (k = 0; k < n - 2 * m1 + 1; k++)
        {
			//get the extents of the first node-------------------
			for (int s = 0; s < 2 * dimension; s += 2)
			{
				rxmbr[s] =    MAXREAL; rxvmbr[s]=MAXREAL;
				rxmbr[s+1] = -MAXREAL; rxvmbr[s+1]=-MAXREAL;
			}
            for (int l = 0; l < m1+k; l++)
            {
				for (int s = 0; s < 2 * dimension; s += 2)
				{
					rxmbr[s] =   min(rxmbr[s],   smu[l].bounces[s]+smu[l].velocity[s]*my_tree->time);
					rxmbr[s+1] = max(rxmbr[s+1], smu[l].bounces[s+1]+smu[l].velocity[s+1]*my_tree->time);
					rxvmbr[s] =   min(rxvmbr[s],   smu[l].velocity[s]);
					rxvmbr[s+1] = max(rxvmbr[s+1], smu[l].velocity[s+1]);
				}
            }
			//add the margin--------------------------------------
			marg += IntgMargin(rxmbr, rxvmbr, 0, my_tree->T);
			//extents of the second node--------------------------
			for (int s = 0; s < 2 * dimension; s += 2)
			{
				rxmbr[s] =    MAXREAL; rxvmbr[s]=MAXREAL;
				rxmbr[s+1] = -MAXREAL; rxvmbr[s+1]=-MAXREAL;
			}
            for (int l= 0 ; l < n; l++)
            {
				for (int s = 0; s < 2 * dimension; s += 2)
				{
					rxmbr[s] =   min(rxmbr[s],   smu[l].bounces[s]+smu[l].velocity[s]*my_tree->time);
					rxmbr[s+1] = max(rxmbr[s+1], smu[l].bounces[s+1]+smu[l].velocity[s+1]*my_tree->time);
					rxvmbr[s] =   min(rxvmbr[s],   smu[l].velocity[s]);
					rxvmbr[s+1] = max(rxvmbr[s+1], smu[l].velocity[s+1]);
				}
            }
			//add the margin--------------------------------------
			marg += IntgMargin(rxmbr, rxvmbr, 0, my_tree->T);
			//----------------------------------------------------
        }

        if (marg < minmarg)
        {
            split_axis = i;
            minmarg = marg;
        }
		//then check velocity dimensions #########################
		//assembly sml and smu for sorting------------------------
        for (int j = 0; j < n; j++)
        {
            sml[j].index=smu[j].index=j;
            sml[j].dimension=smu[j].dimension=i;
            sml[j].bounces=smu[j].bounces=_mbr[j];
			sml[j].velocity=smu[j].velocity=_vmbr[j];
			sml[j].st=smu[j].st=my_tree->time;
			sml[j].ed=smu[j].ed=my_tree->time;
        }
        // Sort by lower and upper -------------------------------
      	qsort(sml, n, sizeof(TPRSortMbr), SortHighVbr);
        qsort(smu, n, sizeof(TPRSortMbr), SortLowVbr);
		//--------------------------------------------------------
        marg = 0.0;
        // for all possible distributions of sml------------------
        for (k = 0; k < n - 2 * m1 + 1; k++)
        {
			//get the extents of the first node-------------------
			for (int s = 0; s < 2 * dimension; s += 2)
			{
				rxmbr[s] =    MAXREAL;  rxvmbr[s]=MAXREAL;
				rxmbr[s+1] = -MAXREAL;  rxvmbr[s+1]=-MAXREAL;
			}
            for (int l = 0; l < m1 + k; l++)
            {
				for (s = 0; s < 2*dimension; s += 2)
				{
					rxmbr[s] =   min(rxmbr[s],   sml[l].bounces[s]+sml[l].velocity[s]*my_tree->time);
					rxmbr[s+1] = max(rxmbr[s+1], sml[l].bounces[s+1]+sml[l].velocity[s+1]*my_tree->time);
					rxvmbr[s] =   min(rxvmbr[s],   sml[l].velocity[s]);
					rxvmbr[s+1] = max(rxvmbr[s+1], sml[l].velocity[s+1]);
				}
			}
			//add the margin--------------------------------------
			marg += IntgMargin(rxmbr, rxvmbr, 0, my_tree->T);
			//extents of the second node--------------------------
			for (s = 0; s < 2 * dimension; s += 2)
			{
				rxmbr[s] =    MAXREAL; rxvmbr[s]=MAXREAL;
				rxmbr[s+1] = -MAXREAL; rxvmbr[s+1]=-MAXREAL;
			}
            for (int l = 0; l < n; l++)
            {
				for (s = 0; s < 2 * dimension; s += 2)
				{
					rxmbr[s] =   min(rxmbr[s],   sml[l].bounces[s]+sml[l].velocity[s]*my_tree->time);
					rxmbr[s+1] = max(rxmbr[s+1], sml[l].bounces[s+1]+sml[l].velocity[s+1]*my_tree->time);
					rxvmbr[s] =   min(rxvmbr[s],   sml[l].velocity[s]);
					rxvmbr[s+1] = max(rxvmbr[s+1], sml[l].velocity[s+1]);
				}
            }
			//add the margin--------------------------------------
			marg += IntgMargin(rxmbr, rxvmbr, 0, my_tree->T);
			//----------------------------------------------------
        }

        // for all possible distributions of smu------------------
       	for (k = 0; k < n - 2 * m1 + 1; k++)
        {
			//get the extents of the first node-------------------
			for (int s = 0; s < 2 * dimension; s += 2)
			{
				rxmbr[s] =    MAXREAL; rxvmbr[s]=MAXREAL;
				rxmbr[s+1] = -MAXREAL; rxvmbr[s+1]=-MAXREAL;
			}
            for (int l = 0; l < m1+k; l++)
            {
				for (s = 0; s < 2 * dimension; s += 2)
				{
					rxmbr[s] =   min(rxmbr[s],   smu[l].bounces[s]+smu[l].velocity[s]*my_tree->time);
					rxmbr[s+1] = max(rxmbr[s+1], smu[l].bounces[s+1]+smu[l].velocity[s+1]*my_tree->time);
					rxvmbr[s] =   min(rxvmbr[s],   smu[l].velocity[s]);
					rxvmbr[s+1] = max(rxvmbr[s+1], smu[l].velocity[s+1]);
				}
            }
			//add the margin--------------------------------------
			marg += IntgMargin(rxmbr, rxvmbr, 0, my_tree->T);
			//extents of the second node--------------------------
			for (s = 0; s < 2 * dimension; s += 2)
			{
				rxmbr[s] =    MAXREAL; rxvmbr[s]=MAXREAL;
				rxmbr[s+1] = -MAXREAL; rxvmbr[s+1]=-MAXREAL;
			}
            for (l=0 ; l < n; l++)
            {
				for (s = 0; s < 2 * dimension; s += 2)
				{
					rxmbr[s] =   min(rxmbr[s],   smu[l].bounces[s]+smu[l].velocity[s]*my_tree->time);
					rxmbr[s+1] = max(rxmbr[s+1], smu[l].bounces[s+1]+smu[l].velocity[s+1]*my_tree->time);
					rxvmbr[s] =   min(rxvmbr[s],   smu[l].velocity[s]);
					rxvmbr[s+1] = max(rxvmbr[s+1], smu[l].velocity[s+1]);
				}
            }
			//add the margin--------------------------------------
			marg += IntgMargin(rxmbr, rxvmbr, 0, my_tree->T);
			//----------------------------------------------------
        }

        if (marg < minmarg)
        {
            split_axis = i+dimension;
            minmarg = marg;
        }
    }
	//end of choosing split axis ---------------------------------
	int dist=-1;
	bool vel_split, lu=true; 
	if (split_axis<dimension)
	{
		vel_split=false;
	}
	else
	{
		vel_split=true;
		split_axis-=dimension;
	}
	//------------------------------------------------------------
    for (int j = 0; j < n; j++)
    {
		sml[j].index = smu[j].index = j;
		sml[j].dimension = smu[j].dimension = split_axis;
		sml[j].bounces=smu[j].bounces=_mbr[j];
		sml[j].velocity=smu[j].velocity=_vmbr[j];
		sml[j].st=smu[j].st=my_tree->time;
		sml[j].ed=smu[j].ed=my_tree->time;
    }

    // Sort by lower and upper -----------------------------------
	if (vel_split)
	{
		qsort(sml, n, sizeof(TPRSortMbr), SortHighVbr);
		qsort(smu, n, sizeof(TPRSortMbr), SortLowVbr);
	}
	else
	{
		qsort(sml, n, sizeof(TPRSortMbr), IntgSortHighMbr);
		qsort(smu, n, sizeof(TPRSortMbr), IntgSortLowMbr);
	}

    float minarea = MAXREAL;
    // for all possible distributions ----------------------------
    for (int k = 0; k < n - 2 * m1 + 1; k++)
    {
		//first check sml ----------------------------------------
		float area=0;
		//get the extents of the first node-----------------------
		for (int s = 0; s < 2 * dimension; s += 2)
		{
			rxmbr[s] =    MAXREAL;  rxvmbr[s]=MAXREAL;
			rxmbr[s+1] = -MAXREAL;  rxvmbr[s+1]=-MAXREAL;
		}
        for (int l = 0; l < m1 + k; l++)
        {
			for (s = 0; s < 2*dimension; s += 2)
			{
				rxmbr[s] =   min(rxmbr[s],   sml[l].bounces[s]+sml[l].velocity[s]*my_tree->time);
				rxmbr[s+1] = max(rxmbr[s+1], sml[l].bounces[s+1]+sml[l].velocity[s+1]*my_tree->time);
				rxvmbr[s] =   min(rxvmbr[s],   sml[l].velocity[s]);
				rxvmbr[s+1] = max(rxvmbr[s+1], sml[l].velocity[s+1]);
			}
		}
		//get the area -------------------------------------------
		area+=IntgArea(rxmbr, rxvmbr, 0, my_tree->T);
		//extents of the second node--------------------------
		for (s = 0; s < 2 * dimension; s += 2)
		{
			rxmbr[s] =    MAXREAL; rxvmbr[s]=MAXREAL;
			rxmbr[s+1] = -MAXREAL; rxvmbr[s+1]=-MAXREAL;
		}
        for (; l < n; l++)
        {
			for (s = 0; s < 2 * dimension; s += 2)
			{
				rxmbr[s] =   min(rxmbr[s],   sml[l].bounces[s]+sml[l].velocity[s]*my_tree->time);
				rxmbr[s+1] = max(rxmbr[s+1], sml[l].bounces[s+1]+sml[l].velocity[s+1]*my_tree->time);
				rxvmbr[s] =   min(rxvmbr[s],   sml[l].velocity[s]);
				rxvmbr[s+1] = max(rxvmbr[s+1], sml[l].velocity[s+1]);
			}
        }
		//get the area -------------------------------------------
		area+=IntgArea(rxmbr, rxvmbr, 0, my_tree->T);
		//check if this distribution is better -------------------
		if (area<minarea)
        {
            minarea = area;
			dist = m1+k;
            lu = true;
        }


		//Then check smu -----------------------------------------
		area=0;
		//get the extents of the first node-----------------------
		for (s = 0; s < 2 * dimension; s += 2)
		{
			rxmbr[s] =    MAXREAL;  rxvmbr[s]=MAXREAL;
			rxmbr[s+1] = -MAXREAL;  rxvmbr[s+1]=-MAXREAL;
		}
        for (l = 0; l < m1 + k; l++)
        {
			for (s = 0; s < 2*dimension; s += 2)
			{
				rxmbr[s] =   min(rxmbr[s],   smu[l].bounces[s]+smu[l].velocity[s]*my_tree->time);
				rxmbr[s+1] = max(rxmbr[s+1], smu[l].bounces[s+1]+smu[l].velocity[s+1]*my_tree->time);
				rxvmbr[s] =   min(rxvmbr[s],   smu[l].velocity[s]);
				rxvmbr[s+1] = max(rxvmbr[s+1], smu[l].velocity[s+1]);
			}
		}
		//get the area -------------------------------------------
		area+=IntgArea(rxmbr, rxvmbr, 0, my_tree->T);
		//extents of the second node--------------------------
		for (s = 0; s < 2 * dimension; s += 2)
		{
			rxmbr[s] =    MAXREAL; rxvmbr[s]=MAXREAL;
			rxmbr[s+1] = -MAXREAL; rxvmbr[s+1]=-MAXREAL;
		}
        for (; l < n; l++)
        {
			for (s = 0; s < 2 * dimension; s += 2)
			{
				rxmbr[s] =   min(rxmbr[s],   smu[l].bounces[s]+smu[l].velocity[s]*my_tree->time);
				rxmbr[s+1] = max(rxmbr[s+1], smu[l].bounces[s+1]+smu[l].velocity[s+1]*my_tree->time);
				rxvmbr[s] =   min(rxvmbr[s],   smu[l].velocity[s]);
				rxvmbr[s+1] = max(rxvmbr[s+1], smu[l].velocity[s+1]);
			}
        }
		//get the area -------------------------------------------
		area+=IntgArea(rxmbr, rxvmbr, 0, my_tree->T);
		//check if this distribution is better -------------------
		if (area<minarea)
        {
            minarea = area;
			dist = m1+k;
            lu = false;
        }
    }

	// take the best distribution --------------------------------
    *_distribution = new int[n];
    for (i = 0; i < n; i++)
    {
        if (lu)
            (*_distribution)[i] = sml[i].index;
        else
            (*_distribution)[i] = smu[i].index;
    }
	//------------------------------------------------------------

    delete [] sml;
    delete [] smu;
    delete [] rxmbr;
//  delete [] rymbr;
	delete [] rxvmbr;
//	delete [] ryvmbr;

    return dist;
}
//*/
/*****************************************************************
This function splits the current node into itself and sn
  
para:
sn: the pointer to the new node 

*****************************************************************/

void RTNode::split(RTNode *sn)
{
	int i;
    int n = num_entries;

    float **mbr_array = new floatptr[n];
	float **vmbr_array = new floatptr[n];
    for (int i = 0; i < n; i++)
	{
       	mbr_array[i]=entries[i].bounces;
		vmbr_array[i]=entries[i].velocity;
	}

	int *distribution;
    int dist = split(mbr_array, vmbr_array, &distribution);

	Entry *new_entries1, *new_entries2;
    new_entries1 = new Entry[capacity];
    new_entries2 = new Entry[capacity];

	for (i = 0; i < capacity; i ++)
	{
		new_entries1[i].init_entry(dimension, my_tree);
		new_entries2[i].init_entry(dimension, my_tree);
	}

    for (i = 0; i < dist; i++)
       	new_entries1[i] = entries[distribution[i]];

    for (i = dist; i < n; i++)
       	new_entries2[i-dist] = entries[distribution[i]];

    for (i = 0; i < n; i++)
    {
       	entries[i].son_ptr = NULL;
       	sn->entries[i].son_ptr = NULL;
    }
    delete [] entries;
    delete [] sn->entries;

    entries = new_entries1;
    sn->entries = new_entries2;

    num_entries = dist;
    sn->num_entries = n - dist;

    delete []mbr_array;
	delete []vmbr_array;
	delete []distribution;
}
//------------------------------------------------------------
void RTNode::write_to_buffer(char *buffer)
{
    int i, j, s;

    // Level
    memcpy(buffer, &level, sizeof(char));
    j = sizeof(char);

    // num_entries
    memcpy(&buffer[j], &num_entries, sizeof(int));
    j += sizeof(int);

    s = entries[0].get_size();
    for (i = 0; i < num_entries; i++)
    {
    	entries[i].write_to_buffer(&buffer[j]);
       	j += s;
    }
}

/*****************************************************************
this function gets the mbr of the entry at some future time, assuming
the information is written at reference time 0

para:
mbr: the mbr at time 0
vbr: the vbr
time: future time at which the mbr is computed
dimension: the dimensionality
*****************************************************************/
void future_mbr(float *_mbr, float *_vbr, float _time, int _dim)
{
	for (int i=0; i<2*_dim; i++)
	{
		_mbr[i]=_mbr[i]+_vbr[i]*_time;
	}
}


/*****************************************************************
this function calculates the integrated area of a moving rectangle
in future time [0, H].
#####NOTE#####
THE CURRENT IMPLEMENTATION WORKS FOR 2D ONLY
##############

para:
mbr: the extents of the rectangle
vmbr: the velocities of the rectangle
st: the starting timestamp
ed: the ending timestamp

*****************************************************************/

float IntgArea(float *_mbr1, float *_vmbr, float _st, float _ed)
{
	float mbr[4];
	for (int i=0; i<4; i++)
	{
		mbr[i]=_mbr1[i]+_vmbr[i]*_st;
	}

	float _H=_ed-_st;

	float s1=mbr[1]-mbr[0];
	float s2=mbr[3]-mbr[2];
	float v1=_vmbr[1]-_vmbr[0];
	float v2=_vmbr[3]-_vmbr[2];

	float ret=s1*s2*_H;
	ret+=(s1*v2+s2*v1)*pow(_H, 2)/2;
	ret+=v1*v2+pow(_H, 3)/3;
	return ret;
}

/*****************************************************************
this function calculates the integrated margin of a moving rectangle
in future time [0, H].
#####NOTE#####
THE CURRENT IMPLEMENTATION WORKS FOR 2D ONLY
##############

para:
mbr: the extents of the rectangle
vmbr: the velocities of the rectangle
st: the starting timestamp
ed: the ending timestamp

*****************************************************************/

float IntgMargin(float *_mbr1, float *_vmbr, float _st, float _ed)
{
	float mbr[4];
	for (int i=0; i<4; i++)
	{
		mbr[i]=_mbr1[i]+_vmbr[i]*_st;
	}

	float _H=_ed-_st;

	float s1=mbr[1]-mbr[0];
	float s2=mbr[3]-mbr[2];
	float v1=_vmbr[1]-_vmbr[0];
	float v2=_vmbr[3]-_vmbr[2];

	float ret=(s1+s2)*_H;
	ret+=(v1+v2)*pow(_H, 2)/2;
	return ret;
}

/*****************************************************************
this function calculates the integrated center distance 
in future time [0, H].
#####NOTE#####
THE CURRENT IMPLEMENTATION WORKS FOR 2D ONLY
##############

para:
mbr1: the extents of the 1st rectangle
vmbr1: the velocities of the 1st rectangle
mbr2: the extents of the 2nd rectangle
vmbr2: the velocities of the 2nd rectangle
H: the parameter

*****************************************************************/

float IntgCenterDist(float *_mbr_1, float *_vmbr1, float *_mbr_2, float *_vmbr2, 
					 float _st, float _ed)
{
	int i;
	float _H=_ed-_st;
	float _mbr1[4], _mbr2[4];
	for (int i=0; i<4; i++)
	{
		_mbr1[i]=_mbr_1[i]+_vmbr1[i]*_st;
		_mbr2[i]=_mbr_2[i]+_vmbr2[i]*_st;
	}

	float deltax[2], deltav[2];
	for (i=0; i<2; i++)
	{
		deltax[i]=(_mbr1[2*i]+_mbr1[2*i+1])/2-(_mbr2[2*i]+_mbr2[2*i+1])/2;
		deltav[i]=(_vmbr1[2*i]+_vmbr1[2*i+1])/2-(_vmbr2[2*i]+_vmbr2[2*i+1])/2;
	}
	float a=deltav[0]*deltav[0]+deltav[1]*deltav[1];
	float b=deltax[0]*deltax[0]+deltax[1]*deltax[1];
	float c=2*(deltax[0]*deltav[0]+deltax[1]*deltav[1]);

	if (a==0 && b==0) return 0;
	if (a==0) return _H*sqrt(c);
	if (c==0) return _H*_H*sqrt(a)/2;
	float A,B,C,D,E;
	A=2*a*_H+b;B=sqrt(a*_H*_H+b*_H+c);
	C=2*sqrt(a);D=4*a*c-b*b;E=8*sqrt(a*a*a);
	float items[4];
	items[0]=A*B/(4*a);
	items[1]=log(A/C+B)*D/E;
	items[2]=b*sqrt(c)/(4*a);
	items[3]=log(b/C+sqrt(c))*D/E;
	return(items[0]+items[1]-items[2]-items[3]);
}


/*****************************************************************
this function is used in split function for sorting the tprsortmbr array


*****************************************************************/

int IntgSortHighMbr(const void *d1, const void *d2)
{
    TPRSortMbr *s1, *s2;
    float erg;
    int dimension;
	float pos1,pos2;
    s1 = (TPRSortMbr *) d1;
    s2 = (TPRSortMbr *) d2;
    dimension = s1->dimension;
	pos1=s1->bounces[2*dimension+1]+s1->velocity[2*dimension+1]*s1->st;
	pos2=s2->bounces[2*dimension+1]+s2->velocity[2*dimension+1]*s2->st;
    erg = pos1-pos2;
    if (erg < 0.0)
		return -1;
    else if (erg == 0.0)
			return 0;
		 else 
			return 1;
}


/*****************************************************************
this function is used in split function for sorting the tprsortmbr array

*****************************************************************/

int IntgSortLowMbr(const void *d1, const void *d2)
{
    TPRSortMbr *s1, *s2;
    float erg;
    int dimension;
	float pos1,pos2;
    s1 = (TPRSortMbr *) d1;
    s2 = (TPRSortMbr *) d2;
    dimension = s1->dimension;
	pos1=s1->bounces[2*dimension]+s1->velocity[2*dimension]*s1->st;
	pos2=s2->bounces[2*dimension]+s2->velocity[2*dimension]*s2->st;
    erg = pos1-pos2;
    if (erg < 0.0)
		return -1;
    else if (erg == 0.0)
			return 0;
		 else 
			return 1;
}


/*****************************************************************
this function is used in insert function for sorting the tprsortmbr array

*****************************************************************/

int IntgSortCenMbr(const void *d1, const void *d2)
{
    TPRSortMbr *s1, *s2;
    s1 = (TPRSortMbr *) d1;
    s2 = (TPRSortMbr *) d2;

    float e1=0, e2=0;
	e1=IntgCenterDist(s1->bounces, s1->velocity, s1->centermbr, s1->centervmbr, s1->st, s1->ed);
	e2=IntgCenterDist(s2->bounces, s2->velocity, s2->centermbr, s2->centervmbr, s2->st, s2->ed);
    if (e1 < e2)
		return -1;
    else if (e1 == e2)
			return 0;
		 else 
			return 1;
}

/*****************************************************************
this function is used in split function for sorting the tprsortmbr array


*****************************************************************/

int SortLowVbr(const void *d1, const void *d2)
{
    TPRSortMbr *s1, *s2;
    float erg;
    int dimension;
	float pos1,pos2;
    s1 = (TPRSortMbr *) d1;
    s2 = (TPRSortMbr *) d2;
    dimension = s1->dimension;
	pos1=s1->velocity[2*dimension];
	pos2=s2->velocity[2*dimension];
    erg = pos1-pos2;
    if (erg < 0.0)
		return -1;
    else if (erg == 0.0)
			return 0;
		 else 
			return 1;
}


/*****************************************************************
this function is used in split function for sorting the tprsortmbr array

*****************************************************************/

int SortHighVbr(const void *d1, const void *d2)
{
    TPRSortMbr *s1, *s2;
    float erg;
    int dimension;
	float pos1,pos2;
    s1 = (TPRSortMbr *) d1;
    s2 = (TPRSortMbr *) d2;
    dimension = s1->dimension;
	pos1=s1->velocity[2*dimension+1];
	pos2=s2->velocity[2*dimension+1];
    erg = pos1-pos2;
    if (erg < 0.0)
		return -1;
    else if (erg == 0.0)
			return 0;
		 else 
			return 1;
}

/*****************************************************************
this function is gets the number of entries in the subtree rooted
at this node
*****************************************************************/

float RTNode::get_subtree_num_entries()
{
	int num=0;
	if (level==0)
	{
		num=num_entries;
	}
	else
	{
		for (int i=0; i<num_entries; i++)
		{
			RTNode *succ=entries[i].get_son();
			num+=succ->get_subtree_num_entries();
			entries[i].del_son();
		}
	}
	return num;
}

/*****************************************************************
this function updates the MBR and VBR according to its child node
para:
pos: the position of the parent entry in the current node

*****************************************************************/

void RTNode::update_parent_entry(int _pos)
{
	RTNode *succ=entries[_pos].son_ptr;
	float *mbr=succ->get_mbr();
	//float *vbr=succ->get_vmbr(); //0906
	float *vbr=succ->get_vmbr(mbr);
	past_mbr(mbr, vbr, my_tree->time, dimension);
	memcpy(entries[_pos].bounces, mbr, sizeof(float)*2*dimension);
	memcpy(entries[_pos].velocity, vbr, sizeof(float)*2*dimension);
	delete [] mbr; delete []vbr;
	dirty=true;
}

/*****************************************************************
this function deletes an entry from the current tree
para:
olde: the entry to be inserted -- need its mbr at the current time set
*****************************************************************/

R_DELETE RTNode::delete_entry(Entry *_olde)
{
	int j;
	R_DELETE ret=NOTFOUND;
	if (level==0)
	{
		for (int i=0; i<num_entries; i++)
		{
			if (_olde->son==entries[i].son)
			{
				dirty=true;
				for (int j=i; j<num_entries-1; j++)
					entries[j]=entries[j+1];
				num_entries--;
				if (num_entries<(capacity-1)*0.4)
				{
					for (j=0; j<num_entries; j++)
					{
						future_mbr(entries[j].bounces, entries[j].velocity,
							my_tree->time, dimension);
						Linkable *link=entries[j].gen_Linkable(); 
						link->level=level;
						my_tree->deletelist->insert(link);
					}
					ret=ERASED;
					my_tree->num_of_dnodes--;
				}
				else
					ret=NORMAL;
				break;
			}
		}
	}
	else
	{
		for (int i=0; i<num_entries; i++)
		{
			float mbr[4];
			memcpy(mbr, entries[i].bounces, 4*sizeof(float));
			future_mbr(mbr, entries[i].velocity, my_tree->time, 2);
			float *ovr=NULL, *vovr=NULL;
			ovr=overlapRect(2, mbr, _olde->bounces);
			//vovr=overlapRect(2, entries[i].velocity, _olde->velocity); //checking overlap velocity, just for traditional MBR
			//if (ovr && vovr)
			if (ovr) //0907 with new bounding rectangle, it is wrong if checking vovr
			{
				delete [] ovr;
				RTNode *succ=entries[i].get_son();
				R_DELETE cret=succ->delete_entry(_olde);
				if (cret==NORMAL)
					update_parent_entry(i);
				entries[i].del_son();
				if (cret==NORMAL)
				{
					ret=NORMAL; break;
				}
				else if (cret==ERASED)
				{
					dirty=true;
					for (int j=i; j<num_entries-1; j++)
						entries[j]=entries[j+1];
					num_entries--;
					if ((level<my_tree->root_ptr->level && num_entries<(capacity-1)*0.4) ||
						num_entries==1)
					{
						for (j=0; j<num_entries; j++)
						{
							future_mbr(entries[j].bounces, entries[j].velocity,
								my_tree->time, dimension);
							Linkable *link=entries[j].gen_Linkable(); 
							link->level=level;
							my_tree->deletelist->insert(link);
						}
						ret=ERASED;
						my_tree->num_of_inodes--;
//printf("underflow in non-leaf level\n");
					}
					else
						ret=NORMAL;
					break;
				}
			}
			else
			{
				delete []ovr; delete []vovr;
			}
		}
	}
	return ret;
}

